Astronomy & Astrophysics manuscript no. mueller201 1 


©ESO 2012 


February 23, 2012 





Circumstellar disks in binary star systems 

Models for y Cephei and a Centauri 

Tobias W. A. Miiller and Wilhelm Kley 



Institut fur Astronomie & Astrophysik, Universitat Tubingen, Auf der Morgenstelle 10, 72076 Tubingen, Germany 
e-mail: Tobias_Mueller@twam. info 



(N 

o 

(N 
X> 

PL, 
(N 



Oh 

43 
6 

> 

in 

oo 



X 



Received 4 October 2011 / Accepted 2 December 201 1 



ABSTRACT 



Context. As of today, over 50 planetary systems have been discovered in binary star systems, some of which have binary separations 
that are smaller than 20 AU. In these systems the gravitational forces from the binary have a strong influence on the evolution of the 
protoplanetary disk and hence the planet formation process. 

Aims. We study the evolution of viscous and radiative circumstellar disks under the influence of a companion star. We focus on the 
eccentric y Cephei and a Centauri system as examples and compare disk quantities such as disk eccentricity and precession rate to 
previous isothermal simulations. 

Methods. We performed two-dimensional hydrodynamical simulations of the binary star systems under the assumption of coplanarity 
of the disk, host star and binary companion. We used the grid-based, staggered mesh code FARGO with an additional energy equation 
to which we added radiative cooling based on opacity tables. 

Results. The eccentric binary companion perturbs the disk around the primary star periodically. Upon passing periastron, spirals arms 
are induced that wind from the outer disk towards the star. In isothermal simulations this results in disk eccentricities up to e iiA as 0.2, 
but in more realistic radiative models we obtain much smaller eccentricities of about ss 0.04 - 0.06 with no real precession. 
Models with varying viscosity and disk mass indicate that disks with less mass have lower temperatures and higher disk eccentricity. 
Conclusions. The fairly high disk eccentricities, as indicated in previous isothermal disk simulations, implied a more difficult planet 
formation in the y Cephei system caused by the enhanced collision velocities of planetesimals. We have shown that under more 
realistic conditions with radiative cooling the disk becomes less eccentric and thus planet formation may be made easier. However, we 
estimate that the viscosity in the disk has to very low, with a < 0.001, because otherwise the disk's lifetime will be too short to allow 
planet formation to occur along the core instability scenario. We estimate that the periodic heating of the disk in eccentric binaries 
will be observable in the mid-IR regime. 
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1. Introduction 



protoplanetary disks - hydrodynamics - radiative transfer - methods: numerical - planets 



Presently, about 50 planets are known to reside in binary stars 
systems. In all systems with solar-type stars the planet orbits 
one of the stars while the secondary star acts as a perturber , 
i.e. these are in a so-called S-type configuration (Dvorak 1986). 
Recently, planets seem to have been detected in evolved binary 
star systems that are in a c ircumbinary (P-type) configuration 
(e.g. lBeuermann et al.ll201 lh . The main observational character- 
istic of the known planets in binary sta r s system have be e n sum - 
marized by Eggenberger et al.l (120041) : iRaghavan et al.l (120061) : 
iDesidera & Barbieril (l2007h . Most planets are in binaries with 
very large separations with semi-major axis beyond 100 AU, 
in particular when detected by direct imaging. Observationally, 
there is evidence for fewer pl anets in binaries with a separa- 
tion of less than about 100 AU dEggenberger et al1l2007l) . in ac- 
cordance with the expectation that binarity constitutes a chal- 
lenge to the planet formation process. Despite this, there are 
several sy stems with a quite c lose binary separation such as 
y Cephei dCampbell et ailll988h. Gliese 86 dOueloz et al.ll2000t 



lEls et al.l200ll). HP 4 1004 dZucker et al1l2004 . and HP 196885 



dCorreia et al 



20081) 



As known from several theoretical studies, the presence of 
the secondary renders the planet formation process more difficult 



than around single stars. Owing to the gravitational disturbance, 
the protoplanetary disk is hotter and dynamically more excited 
such that the coagulation and growth process of planetesimals 
as well as the gravitational instabi lity process in the early phase 
of planet fo rmation are hindered (Nelson l2000t iThebault et al.l 
2004, 2006). This is particularly true for the mentioned closer bi- 
naries with an orbital separation of about 20 AU. Hence this the- 
oretical challenge has put these tighter bin aries into the focus of 
studies on planet formation in binary stars ([Quintana et al.ll2007r, 
iPaardekoooer et al.ll2008"l iKlev & Nelsonll2008l) and on "the ir dy- 
namical stability (IPvorak et al.ll2004T lHaghighipour 2006). The 
most famous example in this category is the y Cephei system, 
which ha s been known for over 20 years to con tain a pro- 
toplanet dCampbell et al.l 119881: lHatzes et all 120031) . In recent 
years the orbital parameter of this s ystem have been updated 
dNeuhauser et al.l 120071 iTorre si 120071) and the basic binary pa- 
rameters are quoted in Table [2] (below) because y Cephei is our 
standard model. At the end of the paper we will briefly discuss 
the a Centauri system as well. 

Because planet formation occurs in disks, their dynamical 
structure is of crucial importance in estimating the efficiency 
of planetary growth processes. The prime effect of the sec- 
ondary star is the truncation of the disk owing to tidal torques 
dArtvmowicz & Lubowll9 94). The truncation radius depends on 
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the mass ratio of the binary, its eccentricity and the viscosity in 
the disk. An additional effect, namely the excitation of eccentric 
modes in the disk, which has been noticed previously in clos e 
binary stars dWhitehursjll988t lLubowlll99U iKlev etai1l2008h . 
has recently drawn attention in studies of planet harboring bin ary 
stars as well dPaardekooper et al.l2008tlKlev & Nelsonl2008l) . In 
these studies it was shown that despite the high eccentricity of a 
y Cephei type binary the disk became eccentric with an average 
ecce ntricity of about edisk ~ 0.12 and a coherent disk preces- 
sion (Kiev & Nelson 20 081). However, this leav es the question 
of numerical effects dPaardekooper et all 120081) . Subsequently 
the influence of self- gravity of the disk has been analyzed by 
iMarzari et al ] d2009ah . who concluded that the inclusion of self- 
gravity leads to disks with lower eccentricity on average. 

Interestingly, there may be observational evidence for tidal 
interactions between a companion on an eccentric orbit and a cir - 
cumstellar disk around the primary star dvan Boekel et al 
as inferred from variable brightness on longer timescales. The 
possibility to detect an eccentric disk th rough spectral line ob- 
servations has been analyzed recently bv lRegalv et al.l (1201 ll) . 

In all simulations mentioned above, the disk has typically 
been modeled using a fixed radial temperature distribution. This 
simplifies the numerics such that no energy equation has to be 
solved but of course it lacks physical reality. In this work we 
will extend previous studies and consider a much improved treat- 
ment of the energy equation. We will do so within the two- 
dimen siona l, flat-disk app r oxima tion following iD'Angelo et al.l 
(l2003h and IKlev & Cridal d2008l) . Here, the viscous heating of 
the disk and the radiative losses are taken into account. We will 
study the structure and dynamical evolution of disks with differ- 
ent masses and binary parameter, and will analyze the influence 
of numerical aspects, in particular boundary conditions. Our first 
target of interest will be the well-studied system y Cephei and we 
will present results on the a Centauri system subsequently. 

2. Model setup 

We assumed that the complete system of primary star, circum- 
stellar disk and secondary star is coplanar. Consequently, in our 
calculations we assumed a flat two-dimensional disk orbiting the 
primary. The disk is modeled hydrodynamically, under the as- 
sumption that the action of the turbulence can be described by a 
standard viscous stress tensor. 



2.1. Physics and equations 

We used cylindrical coordinates (r, ip, z) centered on the primary 
star where the disk lies in the equatorial, z — plane. Because 
our model is two-dimensional (r, tp), we solved the vertically in- 
tegrated versions of the hydrodynamical equations. In this ap- 
proximation the continuity equation is 

' V-(Sv) = 0, (1) 



dt 



where v = (v,., v^) = (v, Qr) is the velocity, and £ = f ^ p dz 
the surface density. As indicated, in the following we will also 
use v and rQ. for the radial and orbital velocity, respectively. The 
vertically integrated equation of radial motion is then 



at Or Or 

and for the azimuthal component 



«9(5> 2 Q) 



dp ff¥ 



(2) 



(3) 



Here p is the vertically integrated pressure, Y the gravitational 
potential of both stars, and /, and de scribe the rad ial and az- 
imuthal forces due to the disk viscosity dMassetl20 02j). Owing to 
the motion of the primary star around the center of mass of the 
binary, the coordinate system is non-inertial, and indirect terms 
were included in the equations of motion to account for this. 
These are included in the potential m. The gravitational influ- 
ence of the disk on the binary is neglected, and the disk in non- 
self-gravitating. The vertically integrated energy equation reads 



de 

— + V(ev) = - P V ■ v + Q + - Q- 
dt 



(4) 



where e is the internal energy density, Q + the heating source term 
and Q the cooling source term. To obtain a fully determined 
system, we additionally used the ideal gas law 



p = <RYJ = (y- l)e , 



(5) 



where T is the temperature in the midplane of the disk, y the 
adiabatic index and % the universal gas constant divided by 
the mean molecular mass, which can be calculated by fi = 
ks/(pm u ), where is the Boltzmann constant, p the mean 
molecular weight and m u the unified atomic mass unit. 

The adiabatic sound speed c s within the disk is then given as 



(6) 



= Vrc s ,i 



where c Sj j S0 = yp/X is the isothermal sound speed. The vertical 
pressure scale height H is then 



H 



c s 



c s 



hr , 



(7) 



where Qk denotes the Keplerian angular velocity around the pri- 
mary and h the aspect ratio. 

For the heating term Q+ we assumed that this is solely given 
by viscous dissipation, and it then is given by 

G+ = 2^ K + 2 < + + ¥ (VV)2 ' (8) 

where v is the kinematic viscosity and cr denotes the viscous 
stress tensor, to be wr itten in polar coordinates. The viscosity v 
is given by v = ac s H dShakura & Sunvaevlll973l) . 

The cooling term Q describes the radiative losses from the 
lower and upper disk surface, which can be written as 

e_=2cr R — , (9) 

Teff 

where <x R is the Stefan-Boltzmann consta nt and r e ff an effective 
optical depth. We follo wed the approac h o flKlev & Cridal (l2008h 
and write, according to lHubenvl d!990l) . 



Teff 



3 

-T + 



V3 



1 



4T + T n 



(10) 



The optical depth follows from r = J p<dz, that can be approxi- 
mated by r w putt where p and k(jd, T) are evaluated at the disk's 
midplane. The vertical density profile of a disk is approximately 
given by a Gaussian and hence 



^pH. 



(ID 



To account for the drop of opacity with vertical height we intro- 
duced a correction factor c\ and write finally 



(12) 
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# 


regime 


Kn Tern 2 / si 


CI 


b 


/ min 1 Y\. \ 

' nun l J 


/ man 1 V\. 1 

111JX L J 


1 

1 


Ice grains 


? v 1 0~ 4 


u 




u 


1 70 
1 /U 


2 


Sublimation of ice grains 


2 x 10 16 





-7 


170 


210 


3 


Dust grains 


5 x ifr 3 





1 


210 


4.6 x 10 3 p"n 


4 


Sublimation of dust grains 


2 x 10 34 




-9 


4.6 x 10 3 p n 


3000 


5 


Molecules 


2 x 10~ 8 


-) 

3 


3 


3000 


1.1 x 10 4 p^i 


6 


Hydrogen scattering 


1 x lfr 36 


1 


10 


1.1 x 10 4 pir 


3 x 10 4 pw 


7 


bound-free & free-free 


1.5 x 10 20 


1 


5 
2 


3 x 10 4 p* 





Table 1. Details of the vari ous opacity regimes by ty pe, showing the transition temperature and the constants kq, a and b. All values 
are quoted in cgs units. See lLin & Papaloizo u ( 1985) for more details. 



The constant c\ — j is obtained by comparing two-dimensional 
disk mod els with the fully th ree-dimensional calculations as pre- 
sented in iKlev et alj d2009l) . For the Rosseland mean opacity k 
we adopted p ower-law dependencies on tem perature and densit y 
described by iLin & Papaloizoul (1 1985b and iBell & Lir] (11994 . 
where 



(13) 



for various opacity regimes. Each opacity regime is described by 
a minimum temperature r m ; n and maximum temperature r max , 
which depends on the density p. Table [TJ li sts the c onstants ko, 
a and b for each regime of the ILin & Papaloizoul model. The 
temperature and density are taken from the midplane, where the 
density p is obtained from Eq. ( fTTT ). 

In our model we did not consider presently any radiation 
transport within the disk plane. This contribution is potentially 
important when strong gradients in temperature and density oc- 
cur. In our situation this may be the case around the periastron 
phase of the binary. However, in our simulations the observed 
contrast did not seem strong enough and we did not expect a 
large impact on the evolution. In subsequent studies we plan to 
investigate this question further. 

Because we are interested in the global evolution of the disk, 
we measured the disk eccentricity edisk and the disk periastron 
Wdisk by first calculating for each grid cell the eccentricity vector 
e, which is defined by 



VXJ 

CM 



(14) 



where / = r X v is the specific angular momentum, G the grav- 
itational constant, M the total mass and r the relative vector. In 
our two-dimensional case the specific angular momentum only 
has a component in z direction and therefore the eccentricity e 
and the longitude of periastron m follow as 



e = \e\ = + e 2 y 
m — atan2(e v , e x ) . 



(15) 
(16) 



The global disk eccentricity edisk and disk periastron nidisk is then 
calculated by a mass-weighted average over the whole disk, 



disk — 



ETdisk - 



I Herd(pdr 
I 'Lm rdtpdr 



Jr min JO 



"Lrdipdr 



Sr dipdr 



(17) 



(18) 



where the integrals are evaluated by summation over all grid 
cells. 



Primary star mass (M primary ) 


1.4 Mq 


Secondary star mass (M seC ondary) 


0.4 Mq 


Binary semi-major axis (a) 


20 AU 


Binary eccentricity (ebin) 


0.4 


Binary orbital period {P^) 


66.6637 a 


Disk mass (M^) 


0.01 Mq 


Viscosity (a) 


0.01 


Adiabatic index (7) 


7/5 


Mean-molecular weight (p) 


2.35 


Initial density profile (S) 


cc r _1 


Initial temperature profile (T) 


CO- 1 


Initial disk aspect ratio (H/r) 


0.05 


Grid (N r x N v ) 


256 x 574 


Computational domain (R m i n - /? max ) 


0.5- 8 AU 



Table 2. Parameters of the standard model. The top entries re- 
fer to the fixed binary parameter. The following give the disk 
properties, and then the initial disk setup and the computational 
parameter are given. 



2.2. Numerical considerations 



The s imulations were performe d using the FARGO code dMassetl 
l2000h with modifications from iBaruteaul (120081) . The numeri- 
cal method used in FARGO is a staggered mesh finite difference 
method. It uses operator splitting and a first-order integrator to 
update to velocities with the source terms (potential and pressure 
gradients, viscous and centrifugal accelerations). The advective 
terms are treated by a second-order upwind algorithm (Ivan Leeri 
119771) . T o speed up cal culations the code uses the FARGO al- 
gorithm dMassedl2000l) . Th e algorithmic d etails of the FARGO 
code have been described in lMassetl d2000l) . Because the FARGO 
code is based on ZEUS - 2D, the basic techniques are described in 
IStone&Normanl(ll992l) . 

The position of the secondary star is calculated by a fifth- 
order Runge-Kutta algorithm. To smooth shocks we u sed the 
artificial viscosity described by IStone & Norma nl <fl992h in our 
simulations. 

To test the code we checked that we obtained the same results 
when using a corotating coordinate system. In addi tion, several 
test calculations were made using the RH2D code dKley||1989l 
Il999h to assure that we implemented everything correctly. 

To avoid numerical problems, we implemented a surface 
density floor of Zg oor = 10~ 7 x So where So = 2(1 AU)| f= o and a 
temperature floor of rfl oor = 3 K, which is about the temperature 
of the cosmic background radiation. In addition to the physical 
motivation for the floor, very low temperatures only occur in the 
very outer parts of the disk (see e.g. Fig. [3} where there is only 
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very little mass that influences the dynamics of the disk. For the 
same reason, these cells do not contribute to the mass weighted 
average of the disk eccentricity. We performed test simulations 
with a very low temperature floor and obtained identical results. 

For the simulations we typically used two boundary condi- 
tions. The open boundary condition allows outflow, but no inflow 
of material. It is implemented as a zero-gradient outflow condi- 
tion, which reads at /? m ; n as 

Tqj = Ei; 
eoj = ex, 

v 0j - v ij — if v 2j > and V2j, otherwise, 

where the index denotes the first inner radial ghost cells, e.g. 
Siy is the first active cell in radial and the j-th cell in azimuthal 
direction. Alternatively, the reflecting boundary condition denies 
exchange of quantities with the system as is implemented at R mm 

as 

E y = Ely 

eoy = 

VQj = V 2 j 
Viy = 

and similarly at /? max . For the outer radial boundary we always 
used the standard outflow condition, while at /? m ; n we typically 
used the reflecting condition. Below, we will investigate a possi- 
ble influence of the numerical inner boundary condition. 

3. The standard model 

For our simulations we constructed a standard reference model 
using the physical parameter of the y Cephei system. The first 
entries of Table[2]list the orbi tal parameter of this stan dard model 
based on the observations of iNeuhauser et all d2007l) and lTorrell 
d2007h . 

The disk parameters where chosen to be in a typical range 
for a possible disk around y Cephei A, with a total mass that 
would allo w for the format ion of the observed planet of about 
1.85Mj up dEndl et alJl201lh in the system. Because the tidal 
forces of the secondary limits the radial extent of the disk to 
about 6AU dArtvmowicz & Lubowlll994T) we chose a slightly 
larger extent of the computational domain to 8 AU. This reduces 
artificial boundary effects. The initial density distribution is lim- 
ited to 6 AU. For the outer boundary condition we always used 
the zero-gradient outflow boundary condition so that the material 
that is accelerated during periastron may leave the system freely. 
The inner boundary in the standard model is typically reflect- 
ing. The resolution of the logarithmic grid is trimmed to have 
quadratic cells (rAtf/Ar as 1). 

After the initial setup at time zero, the models have to run 
for several binar y orbits until a quasis tationary state has been 
reached, see also lKlev & Nelsonld2008l) . This equilibration pro- 
cess typically takes about 15 binary periods, i.e. 1000 years for 
fully radiative models. Afterward the disk cycles through sev- 
eral states during one orbital period of the binary displayed in 
Fig. [T] When the binary is at apastron, the the disk is very ax- 
isymmetric, but when the binary moves toward its periastron, it 
starts to perturb the disk. After the binary passes its periastron, 
two strong tidal spiral arms develop within the disk and wind 
themselves to the center of the disk. Before the binary reaches 
its apastron, the disk spiral arms disappear. Periastron passages 
are also the time when most of the disk mass is lost. During the 
first periapse the disk looses 4.8 % of its mass. This decreases 



during the next ten periapsis passes to about 0. 1 %. At the end of 
the simulation after 200 binary orbits, the disk mass has reduced 
to 0.0073 Mq, which is a loss of about 20 %. 

3. 1 . Isothermal runs 

To check our results and compare them with previous results, we 
first performed simulations where we kept the initial temperature 
stratification. In Fig. [2] we present results for constant H/r disks 
where we did not solve the energy equation. Shown are results 
for four different values of H/r, which were chosen according 
to our fully radiative results (see below). The disk reaches an 
equilibrium state after many (typically several 100) binary or- 
bits, and the magnitude of the final average eccentricity depends 
on the temperature of the disk. In particular, for the disks with 
an H/r of about 0.05 the disk eccentricity settles at fairly high 
values of about 0.2. Thi s seems to be higher t han in previous 
simulations performed bv lKlev & Nelson! d2008l) . but the param- 
eters used for the y Cephei system, in particular the mass ratio 

q = Msecondary/Afpnmary °f tne nost star an d me binary companion 
(their q = 0.24 instead of our q = 0.29), and the binary eccen- 
tricity ebin (<?bin = 0.36 instead of e^n = 0.4 ) are differe nt. Using 
the (older) system parameters as given by Kiev et afl. we were 
able to reproduce the results nicely (e.g. their Fig. 4). 

In these isothermal runs it is clear that the two thicker, hotter 
disks have significantly lower mean eccentricities. The disk with 
H/r = 0.065 does not even display a coherent disk precession 
any more (Fig. [2] bottom panel). This drop of edisk with increas- 
ing H/r is a consequence of the reduction of the eccentricity 
growth rate for thicker disks. The growth rate of the disk eccen- 
tricity is the lowest for the hottest disk with H/r = 0.065 and 
increases with decrea sing H/r (F i g. [2] t op panel), a finding that 
perfectly agrees with iKlev et all d2008l) . The hotter disks need 
therefore much more time to reach an equilibrium state. As the 
eccentric binary damps the growth of disk eccentricity by tidal 
interaction (see below), the hotter disks settle into an equilibrium 
state that has a lower edisk- 

After about 400 binary periods, i.e. 25000 years, the time av- 
erage of the disk eccentricity (edisk) reaches a constant value and 
the disks precess with a nearly constant precession rate crdisk- 
The values of (edisk) and crdisk for the four different values of 
H/r are given in Table [3] 



H/r 


(edisk) 


ra'disk [P hiB l 


0.05 


0.19803 ±0.00085 


-0.0844 ±0.0018 


0.055 


0.20200 ± 0.00076 


-0.1082 ±0.0017 


0.06 


0.15800 ±0.00054 


-0.1 172 ±0.0063 


0.065 


0.08004 ±0.00031 


-0.1334 ±0.0025 



Table 3. Time average of the disk eccentricity (edisk) and pre- 
cession rate riidisk for the four different isothermal models. 



3.2. Fully radiative disks: Structure and dynamics 

Now we include all physics, in particular the viscous heating 
and radiative cooling terms and solve the full set of equations as 
stated above. For given binary parameter and disk physics (opac- 
ity and viscosity) the dynamical state and structure of the disk is 
completely determined once the disk mass has been specified. 
The density and temperature distribution cannot be stated freely 
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-8 -6 -4 -2 2 4 6 8 -8 -6 -4 -2 2 4 6 8 -8 -6 -4 -2 2 4 6 8 -8 -6 -4 -2 2 4 6 8 



x [AU] x [AU] x [AU] x [AU] 

Fig. 1. Evolution of the disk within one binary period for the fully radiative standard model. The top labels in the four panels state 
the current position of the binary (rbi n , <fbm) in polar coordinates. In the first panel the binary companion is at the apastron and 
has less influence on the disk. The L\ point in this configuration is at about 17.5 AU and therefore far outside the computational 
domain (gray circles). In the second panel the binary has reduced its distance to about 12.8 AU and the Roche lobe (green curve) 
has shrunken dramatically and is now entirely within the computational domain. The third panel shows the disk with the binary at 
periastron. Some of the material in the disk is now outside the Roche lobe and might be lost from the system. In the last panel the 
binary's separation increases again, which makes the Roche lobe grow such that it engulfs the disk entirely. The strong tidal forces 
near periastron induce spiral waves in the disk that will be damped out until the binary is at apastron. 



Time [P bin ] 

50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 




5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 
Time [a] 
Time [P bin ] 

10 20 30 40 50 60 70 80 90 100 




Fig. 2. Global mass-weighted disk eccentricity {top) and disk 
periastron {bottom), both sampled at the binary's apastron, for 
isothermal simulations of the standard model with different as- 
pect ratios H/ r. The disk eccentricity needs quite a long time to 
reach a quasi-equilibrium and settles at around 0.2 for an H/ r of 
0.05 and 0.055. High values of H/r result in a lower disk eccen- 
tricity. Note that in the lower panel a short time span is displayed 
for clarity. 




0.14 
0.12 
0.1 
0.08 
0.06 
0.04 
0.02 




= 0Pbi, 
-- 25 P bjn = 1667 a 
= 50 P bin = 3333 a 
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Distance [AU] 




= 0P„„ 
= 25 P bin = 1667 a 
= 50 P bin = 3333 a 
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Distance [AU] 



t = 0P bi „ = 0a 
t = 25 P bjn = 1667 a 
t = 50 P bin = 3333 a 
t = 75P bm = 5000a 
I = 100 P bi „ = 6666 a 



4 5 6 

Distance [AU] 



Fig. 3. Radial dependency of surface density {top), midplane 
temperature {middle) and disk eccentricity {bottom) of the fully 
radiative standard model for different timestamps. The disk 
ranges from 0.5 AU to about 5.5-6 AU due to tidal forces of the 
companion. At t — years the disk is not fully Keplerian and 
thus e disk > 0. 
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Fig. 4. Radial dependency of surface density (top), and disk ec- 
centricity (bottom) of the H/r = 0.05 isothermal model and the 
fully radiative standard model after 750 binary orbits. 



and hence the power law given in Table|2]refers only to the initial 
setup. 

The radial disk structure is displayed in Fig. [3] where az- 
imuthally averaged quantities of S, T and edisk are displayed. The 
upper panel of Fig.[3]shows the surface density and temperature 
distribution for five different timestamps. For the initial setup, 
at t = years, the surface density follows the initial r profile 
in the inner region (0.5 AU-5 AU) and is lowered to the density 
floor in the outer region. The temperature follows the initial r _1 
profile for the whole computational domain. After about 25 bi- 
nary orbits a quasistationary state is reached where the profiles 
show a bend at about r — 1.8 AU. This change in slope, which 
occurs at a temperature of about 1000K in the profiles, is a re- 
sult of a change in the opacity caused by the sublimation of dust 
grains, which starts at about (pg^'cm 3 ) 1 ^ 15 x 4.6 x 10 3 K in th e 
opacity tables (see Table [TJ given bv lLin & Papaloizoul d!985h . 
A second bend occurs at about r = 4.5 AU, which corresponds 
to the truncation radius of the companion's tidal forces. In the 
bottom panel of Fig. [3] the corresponding radial distribution of 
eccentricity is shown. The eccentricity is low in the inner parts 
of the disk and increases with larger radii. The highest values 
for e(r) occur in the outer region r > 5.5 AU, beyond the tidal 
truncation radius. But these regions do not play an important dy- 
namical role because there is only very little mass left owing to 
the tidal forces of the secondary star. For example, the disk has 
a mass-weighted mean overall disk eccentricity of edisk = 0.032 
after 100 binary orbits even though more than half of the ra- 
dial extent of the disk has an eccentricity of > 0.05. Remarkably 
the eccentricity is much lower than in the isothermal simulations 
presented in Section lXTl As a consequence of the lack of eccen- 
tricity there is not visible disk precession. 

Fig. [3] shows the radial dependency of the disk surface den- 
sity and disk eccentricity of the H/r = 0.05 model in equilibrium 
after 750 binary orbits compared with the subsequently intro- 
duced fully radiative model. The surface density in the isother- 
mal model misses the bend caused by the 1000K change in the 




1 2 
Distance [AU] 

Fig. 5. Azimuthally averaged surface density (top) and midplane 
temperature (bottom) profiles for the radiative standard model 
using two different locations of the inner boundary of the com- 
putational domain. Results are displayed at 100 binary orbits 
(6666 years). In the overlapping region from 0.5-8 AU the pro- 
files match very well and the region from 0.25-0.5 AU in the 
'"min = 0.25 AU model is a consistent continuation. 



disk opacity and is slightly steeper than the initial r profile. 
The disk eccentricity in the isothermal case is much more homo- 
geneous with radius, as is to be exp ected for a coherent str ucture. 
This has already been observed bv lKlev & Nelson! (l2008h . 

To check for numerical dependencies we varied the loca- 
tion of the inner boundary of the computational domain and 
the grid resolution. First we chose a smaller inner boundary of 
'"min = 0.25 AU increasing the number of radial grid points such 
that the spatial grid resolution remained unchanged in the over- 
lap region. In Fig.[5]the radial dependency of the surface density 
and the midplane temperature at the end of the simulations dif- 
fers not too much in the overlap region of both simulations, so 
that the effects of the inner boundary are negligible in this re- 
spect. However, the amplitude of the oscillation in the time evo- 
lution of the disk eccentricity is nearly gone (see Fig. |6). In the 
standard model we could see a small oscillation in the disk pe- 
riastron, which also vanishes in the r m ; n = 0.25 AU model. The 
non-existence of a coherent precession is also an indication of 
very low disk eccentricity in both cases. In additional test simu- 
lations (not displayed here), we extended the outer radius to the 
value of r max = 12 AU, changing the number of grid cells ac- 
cordingly to maintain the resolution of the standard model. As 
expected, this larger radial domain does not change the results. 

Then we varied the grid resolution by first doubling the grid 
cells to 512 x 1150 and doubling it again to 1024 x 2302 grid 
cells. This change in resolution has basically no impact on the 
disk's dynamics. As before, the disk eccentricity settles to a time 
averaged value of about 0.04. In all three simulations the disk 



6 



T. W. A. Miiller and W. Kley: Circumstellar disks in binary star systems 



Time [P bin ] 
50 75 




1000 2000 3000 4000 5000 6000 7000 8000 

Time [a] 
Time [P bin ] 
50 75 




1000 2000 3000 4000 5000 6000 7000 8000 

Time [a] 

Fig. 6. Global mass-weighted disk eccentricity (top) and disk 
periastron (bottom), sampled at the binary's apastron, for two 
different inner boundary positions. 



periastron does not reach a state with real precession and the 
only change caused by the change of resolution is a small shift 
in time. Consequently our standard resolution of 256 x 574 cells 
seems to be sufficient to resolve the disk dynamics properly. 

4. Variation of physical parameters 

To study the influence of the physical conditions on the evolution 
of the disk, we investigated in particular the impact of the disk 
mass Mdi s k, the viscosity v, which is determined by a, the opacity 
k, and the binary's eccentricity ebin- 

4.1. Disk mass 

The first parameter we investigated is the disk mass. In contrast 
to isothermal simulations our radiative simulations depend on 
the disk mass as the opacity depends on gas density. To analyze 
the influence of the disk mass on the disk's evolution, we ran 
four simulations with disk masses M d ; s k of 0.005 Mq, 0.01 Mq, 
0.02 Mq and 0.04 Mq while keeping all other parameters un- 
changed. 

The surface density and temperature profiles after 100 bi- 
nary orbits for different disk masses is displayed in Fig. [7] where 
the solid (red) line refers to the standard model displayed in 
Fig. [3] The profiles can be divided into two regimes that are sep- 
arated by the 1000 K temperature line caused by the opacity (see 
Section [3~2l i. In each regime the surface density and temperature 
profiles follow a simple power-law that depends only weakly on 
the disk mass, and is given in the caption of Fig. [7] 

Fig- E shows the time evolution of the disk eccentricity and 
periastron for different disk masses. The higher the disk mass, 
the lower the oscillations of the the disk eccentricity, but the 



time average of the disk eccentricity is in the range of 0.04 - 
0.05 for all disk masses. The disk periastron displays no real pre- 
cession and with increasing disk mass it settles at about 0. This 
is in contrast to the isothermal simulations presented in Section 
13. II where we obtained a disk eccentricity e d ; s k of 0.2 for an as- 
pect ratio H/r of 0.05 and a real precession of the disk perias- 
tron. There is a trend, however, that the eccentricity becomes 
higher for cooler disks with lower H/ r. But one has to be care- 
ful here, because the aspect ratio is constant in radius and time 
for the isothermal simulations, but not for our radiative simula- 
tions. Fig. [9] shows the aspect ratio of the disk after 100 binary 
orbits. All models had an initial value of (H/r) initi . d i = 0.05 but 
end up with different values of H/ r depending on the disk mass 
and the phase in the binary orbit. We note that in all models the 
mass of the disk reduces with time owing to the mass loss across 
the outer boundary. In particular, we find that after 100 binary 
orbits the disks have lost about 27 % of their initial mass in the 
M disk = 0.04 Mq model, 23 % in the M disk = 0.04 Mq model, 
17 % in the M disk = 0.04 Mq and 6 % in the M disk = 0.005 Mq 
model. Hence, the values quoted in the text and figures always 
refer to the initial disk masses. 

The results for e d isk in Fig. [8] show a marginal increase of 
the disk eccentricity for smaller disk masses. To test if this trend 
continues, we performed additional simulations for even smaller 
initial disk masses and found indeed an increased oscillatory be- 
havior of the eccentricity, which settles eventually after about 
600 binary orbits to a low eccentric state very similar to the 
A^disk = 0.025 model, however. Hence, there does not seem to 
exist an obvious trend of e d ; s k with disk mass. 

4.2. Viscosity 

In this section we investigate the influence of the viscosity v. To 
study the dependence on v we varied the a parameter, which de- 
termi nes the viscosity, v = ac s H — acl O^ 1 ( Sha kura & Sunvaevl 
119731) . from our standard model (a = 0.01)andkept all other pa- 
rameters unchanged. We varied a from 0.005 to 0.04. 

Fig. [10] shows the surface density and temperature profiles 
after 100 binary orbits. All models started with identical disk 
mass, and evidently, higher a models lose the disk masses more 
rapidly. For example, at the displayed time at 100 binary orbits 
the a = 0.04 model has lost about 71 % of its initial mass, while 
the standard model (a = 0.01) lost only 19 %, see also Fig. [T_2] 
below. The shape of the surface density and temperature profile 
seems to be independent of the viscosity in the disk, as expected. 
Again, the profiles can be divided into two regimes that are sep- 
arated by the 1000K temperature line. The power-laws for the 
surface density and temperature profiles are given in the caption 
ofFig.QJJ 

Fig.fTTIshows the influence of a on the disk eccentricity and 
periastron. Higher values of a and thus high viscosities result in 
a calmer disk that does not react fast enough on the disturbances 
of the binary companion. Therefore the disk shows less oscil- 
lations in the disk eccentricity and the disk periastron remains 
almost constant for the higher values of a. The disk eccentricity 
in the model with the highest viscosity, a = 0.04, increases to 
up to 0.25 within 600 binary orbits, and the disk eccentricity in 
the a - 0.02 model starts to grow slowly after about 200 - 300 
binary orbits. This effect of a rising disk eccentricity is a direct 
consequence of the increasing viscosity. As shown in Fig. [12] the 
mass loss of the disk depends on the magnitude of the viscosity, 
the higher a, the faster the mass loss of the disk. For example, 
after 500 binary orbits the a = 0.04 disk has lost about 97 % 
of its initial mass. The increased mass loss for higher viscosities 
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Fig. 7. Azimuthally averaged surface density {top) and mid- 
plane temperature (bottom) profiles after 100 binary orbits 
(6666 years) for radiative models using different initial disk 
masses. The profiles can be fitted using power-laws when di- 
vided into two regimes: An outer region with temperatures of 
less than about 1000 K and an inner region with temperatures 
above 1000 K because of a break in the opacity tables at about 
1000 K. The inner region follows S oc r ~ ■ to ~' 51 and T oc 
r -o.53 to -0.54^ wn ereas the outer region follows S oc r -° &2 t0 - 1 - 67 
and T oc r~ l 39 t0 58 for all models. 



is a result of the stronger outward spreading of the disk, i.e. a 
larger truncation radius. The larger disk makes the outer parts of 
the disk more susceptible to the tidal perturbations of the sec- 
ondary, which increase the disk eccentricity dramatically. This 
is in accord ance with the visc osity dependence found in earlier 
studies, e.g. lKlev et all (120081) . 



4.3. Opacity 

The amount of cooling in our radiative model depends on the 
disk's opacity. To examine the influence of the opacity model 
we calculated the standard model with two different opacity 
models leaving all other parameters unchanged. In the first cal- 
culation we used the op acity tables (see Table Q} shown by 
iLin & Papaloizoul (Il985h . whereas the secon d calculation uses 
the opacity tables given bv lBell & Linl (1 1994 their Table 3). 

As expecte d, the opac ity plays an important role and the 
model with the lBell & Linl opacity shows less oscillations in the 
disk eccentricity, b ut the time averag e is about the same as in 
the model with the ILin & Papaloizoul opacity. This is similar to 



the fmin change in Section [3T2l The surface density and midplane 
temperature profiles at th e end of the simulation in the simula- 
tion with the iBell & Linl opacity model differ slightly from the 
one shown in Fig. [3] Both models have a bend at a tempera- 
ture of about 1000K but the temperature profile in the model 
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Fig. 8. Global mass-weighted disk eccentricity (top) and disk 
periastron (bottom), sampled at the binary's apastron, depen- 
dency for different disk masses. 
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Fig. 9. H/r after 100 binary orbits (6666 years) for different disk 
masses. 



with the lBell & Linl opacity model is flatter and slightly lower in 
the region below 1000 K and slightly steeper in the region about 
1000K. In exchange, the surface density is slightly steeper and 
slightly higher in region the below 1000 K and matches the other 
model in the region above 1 000 K very well. 

4.4. Binary eccentricity 

Another very important factor for the disk's evolution are the 
binary parameters. We therefore varied the binary's eccentricity 
in our standard model from ebin = to 0.4. Because this also 
changes the truncation radius of the disk that is caused by the bi- 
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Fig. 10. Azimuthally averaged surface density {top) and mid- 
plane temperature {bottom) profiles after 100 binary orbits 
(6666 years) using four different values for the viscosity param- 
eter a. The profiles can be fitted using power-laws when di- 
vided into two regimes: An outer region with temperatures of 
less than about 1000 K and an inner region with temperatures 
above 1000 K, because there is a break in the opacity tables at 
about 1000 K. The inner region follows S °c r~ l 48 10 ,53 and 



.-0.53 



whereas the outer region follows S oc r 



-0.78 to -0.95 



and r oc r 1 45 10 1 53 for all a values. 



nary's tidal forces dArtvmowicz & Lubowl 1 199*41) . we extended 
the computational domain to up to 12.5 AU for the e D in = 
model. To reach the same resolution in the computational do- 
main of the standard model (0.5 - 8 AU) we increased the num- 
ber of cells in radial direction to 295 in the e D in — model. The 
<?bin = 0.05, ebin =0.1 and e D i n = 0.2 models were adjusted ac- 
cordingly in their computational domain and resolution in radial 
direction. 

Fig-EUshows the time evolution of the disk eccentricity and 
periastron. Interestingly, the e D in = and e D in = 0.05 models 
show the highest dis k eccentrici t y. The se high values for e D i n = 
seem to agree with iKlev et all (120081) . Also, the e^n = and 
<?bin = 0.05 models are the only ones that have a real coher- 
ent disk precession with a precession rate of -0.033 P bin for 
the ebin = model and -0.062 PjV for the ebin = 0.05 model. 
All other models with e^ n > 0.1 show a disk eccentricity of 
<?disk < 0.1 and no precession, which also indicates that the disk 
is not globally eccentric. For the low eccentric binaries it takes 
up to about 125 binary orbits to reach the high eccentric quasi- 
equilibrium disk state, which is long compared to the standard 
model, which reaches its quasistationary state after only about 
15 binary orbits. The se timescales agree well with those obtained 
bv IKlev etalJd2008h for isothermal disks. The reason why bina- 
ries with low ebin tend to have eccentric disks is the larger disk 
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Fig. 11. Global mass-weighted disk eccentricity {top) and peri- 
astron {bottom), sampled at the binary's apastron, as a function 
of time for different a values for the viscosity. 
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Fig. 12. Disk mass evolution for different a values for the disk 
viscosity. 



radius in this case. This allows an easier oper ation of the insta- 
bility according to the model of lLubowld 19911) . 

The surface density and midplane temperature profiles af- 
ter 100 binary orbits of all five models have the same slope 
but differ slightly in absolute values. The eb; n = 0.4 models is 
the hottest and densest model and then the temperature and sur- 
face density decreases with decreasing binary eccentricity. The 
ebin = model shows low oscillations in the profiles. As ex- 
pected, the disk's truncation radius owing to the binary's tidal 
forces is shifted outward in the models with lower binary eccen- 
tricity. 
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Fig. 13. Global mass-weighted disk eccentricity (top) and disk 
periastron (bottom), sampled at the binary's apastron, for dif- 
ferent binary eccentricities. The low eccentricity (ebin < 0.05) 
models need much more time to reach their equilibrium state 
compared to the standard model. These are also the only models 
that develop a real coherent disk precession. 

5. Brightness variations 

To identify possible observable changes in the brightness of the 
systems caused by perturbations in the disk, we calculated theo- 
retical light curves for our standard model. For that purpose, we 
examined the time variation of the disk dissipation given by 

£>di.sk = jjf Q + dA, (19) 
and the disk luminosity given by 

W = ff Q dA = JJ 2o- R ^ dA . (20) 

To identify the source of the brightness variations we divided the 
disk into five rings ranging from 0.5-2 AU, 2-3 AU, 3-4 AU, 4- 
5 AU and 5-8 AU. In each of these rings the disk dissipation and 
disk luminosity was calculated. Fig.[14]displays the variation of 
the disk dissipation and disk luminosity of our standard model 
during one orbital orbit at t = 100 Pbin where the system has 
already reached its quasistationary state. The periastron occurs 
at t = 100.5 Pbin, and is indicated by the vertical line. 

When the binary is at apastron, the disk is almost uniform 
(see Fig. |T} and we cannot see any variations in the light curves. 
However, when the binary passes periastron, it starts to perturb 



Fig. 14. Variation of disk dissipation and disk luminosity dur- 
ing one binary orbit. The dissipation and luminosity are calcu- 
lated for five different rings ranging from 0.5-2 AU, 2-3 AU, 
3^1 AU, 4-5 AU and 5-8 AU. Because the disk is truncated at 
about 4.5 AU (see Fig. |3J, the outmost ring does not contain 
much mass. The solid (red) curve is the sum of all five rings. 
The gray line indicates the binary's periastron passage. 



the outer regions of the disk and two spiral arms evolve that 
wind themselves to the center of the disk. This results in a rise 
of the luminosity by 2-3 magnitudes in the outer rings of the 
disk shortly after periastron passage (see Fig. [T4l and after a 
short time also in a smaller rise in the inner rings of the disk. 
As the binary moves farther away from the disk, the spirals are 
damped out and the luminosity peak vanishes. These luminos- 
ity peaks should be observable as a 2-magnitude increase in the 
mid-infrared (MIR) because the main contributions come from 
the outer rings (3-5 AU) with temperatures between 150 K and 
450 K. 

Luminosity peaks have alrea dy been observed in th e mid- 
infrared in the T Tau S system (Ivan Boekel et alJ 12010b . a bi- 
nary system that is not very di fferent from the ea rly stage of the 
y Cephei system. Additionally. Ivan Boekel et"al] also performed 
some radiative simulations for the T Tau S system and pointed 
out that these brightness variations could be caused by the pertur- 
bations of the binary companion. However, the brightness varia- 
tions found in their disk models were very weak. 
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Table 4. Parameters of the a Centauri model. The top entries 
refer to the fixed binary parameter. Below we list along with 
the disk properties the initial disk setup and the computational 
parameters. 



Primary star mass (M primary ) 


1.1 Mq 


Secondary star mass (M sccon dar y ) 


0.93 M Q 


Binary semi-major axis (a) 


23.4 AU 


Binary eccentricity (fibin) 


0.52 


Binary orbital period (.Pbin) 


79.4431 a 


Disk mass (M disk ) 


0.01 Mq 


Viscosity (a) 


0.01 


Adiabatic index (7) 


7/5 


Mean-molecular weight (/j) 


2.35 


Initial density profile (£) 


cc r _I 


Initial temperature profile (r) 


cc r -1 


Initial disk aspect ratio (H/r) 


0.05 


Grid (N r x Ay 


256 x 574 


Computational domain (R mm - R mdx ) 


0.5- 8 AU 



6. a Centauri 

In addition to the y Cephei system we also performed some 
calculations for the a Centauri system, because it is of special 
interest, being the nearest star to our solar system. This sys- 
tem has bee n investigated for the possibility of planet forma- 
tion, see e.g. lThebault etd](l2008b . Tableg] gives an overview of 
the pa rameters of our a Centauri model based on lPourbaix et al.l 
(120021) . We tried to keep the disk parameters the same as in the 
y Cephei model, so that the main difference are the mass ratio 
q = Msecondary/A^pi-iinmy and the binary's orbital parameters. In 
the y Cephei model we have a mass ratio q of 0.28, an eccentric- 
ity ^bin of 0.4 and a semi-major axis a of 20 AU, whereas in the 
a Centauri model we have a mass ratio q of 0.84, an eccentricity 
ebin of 0.52 and a semi-major axis a of 23.4 AU. 

In Section l4~4l we saw that for high binary eccentricities the 
disk eccentricity does not reach very high values, so that we 
would e xpect a fai r ly low disk eccentricity for the a Centauri 
system. iKlev et all (120081) showed that for isothermal simula- 
tions the disk eccentricity in the quasistationary state does not 
depend heavily on the mass ratio q. Fig.[T5lshows the evolution 
of the disk eccentricity and periastron over time. The disk eccen- 
tricity settles after about 15 binary orbits at a rather low value of 
about 0.038 and the disk periastron is also at a nearly constant 
position. 

The surface density and midplane temperature profiles for 
the a Centauri model can only be fitted in the inner region with 
r < 1.8 AU as a simple power-law. The surface density can be 
described by a r -1 - 49 and the midplane temperature by a r~ 0,53 
power-law. In the outer region the surface density and tempera- 
ture decreases rather fast to until about 4 AU, where the disk is 
truncated by the binary's tidal forces. 

7. Summary and conclusions 

We have investigated the dynamics of a protostellar disk in bi- 
nary star systems using specifically the orbital parameter of 
y Cephei and a Centauri. We assumed a coplanar system and 
used a two-dimensional hydrodynamical code to evolve the non- 
self-gravitating disk. Extending previous simulations, we in- 
cluded internal viscous heating given by an a-type viscosity pre- 
scription, and radiative cooling from the disk surface. 
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Fig. 15. Global mass-weighted disk eccentricity {top) and disk 
periastron {bottom), sampled at the binary's apastron for the 
a Centauri model. 



In a first set of simulations we investigated locally isother- 
mal disks for different disk temperatures. We showed that disks 
in binaries of the y Cephei type with a standard thickness of 
H/r — 0.05 become eccentric (eaisk ~ 0.2) showing a coherent 
disk precession . This agree s well with previou s simu lations by 
lKlev & Nelsonl (120081) and iPaardekooper et all (120081) . Varying 
the temperature in the disk, we showed that the magnitude of 
the disk's eccentricity becomes lower when the disk thickness 
increases. For disks with H/r > 0.065 the mean average ec- 
centricity has dropped below 0.08 and the disks do not show a 
precession anymore. 

Then we studied more realistic disks with internal heating 
and radiative cooling, varying the disk's mass. In all cases we 
found relatively low eccentricities and no precession. We at- 
tribute the lack of eccentricities firstly to the increased disk 
height, which is, in particular for the more massive disks, higher 
than the standard value (see Fig. |9). Secondly, in the full ra- 
diative models the disk's dynamical behavior is more adiabatic 
compared to the locally isothermal case. Then, through com- 
pressional heating {pdV-woik), kinetic energy is transferred to 
internal energy, which leads to a reduced growth of disk eccen- 
tricity. We have checked that purely adiabatic models show an 
even lower disk eccentricity than the radiative models. Hence, 
the radiative case lies between the adiabatic and isothermal, as 
expected. 

Because the disk's energy balance is determined via the vis- 
cosity, we changed the value of the parameter a ranging from 
0.005 to 0.04, all values that are consistent with the results 
of MHD-turbulent accretion disks. Here, we found that only 
the disk with the highest a becomes eccentric. The reason for 
this rise is the larger disk radius, which leads to an enhance- 
ment of the tidal torques from the secondary. We note that the 
disk's outer radius in our models still lies well inside the 3:1 
resonanc e with the bin ary. According to the linear instability 
model bv lLubowl(ll991l) . the disk eccentricity is excited through 
the 3:1 resonance and hence, the disk should be sufficiently 
large, a condition which is fulfilled only for s mall mass ratios, 
a - M secondarv/Af primary- However, as shown by IKlev & Nelsonl 
(l2008h . disks in binary star systems with large mass ratios can 
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turn eccentric as well, even though the disks are small, a feature 
confirmed in our simulations. 

The inferred short lifetime of disks with standard viscosities 
is slightly alarming with respect to planet formation in these sys- 
tems. In the core-accretion scenario planet formation proceeds 
along a sequence of many steps that take a few Myr. For disks 
to persist this long in y Cephei-type binaries a very low vis- 
cosity of a < 1CT 4 seems to be required. In the gravitational 
instability scenario the timescale for planet formation is much 
shorter and hence, this scenario may be favored by our find- 
ings. Observationally, several recent studies indeed suggest that 
the lifetime of disks in young binary stars i s significantly re- 
duced compared to disks around single stars (ICieza et al.ll2009t 
lDuchend2010tlKraus et al.ll201 ll) . Dynamically, this behavior is 
expected, because the perturbation of the companion star leads 
quite naturally to an increased mass loss of the disk, details de- 
pending on the binary separation and eccentricity. 

A very critical phase in the core accretion scenario, in par- 
ticular in binary stars, is the initial growth of meter to km- 
sized planetesimals. Here, the growth depends on the successful 
sticking of the two collision partners. Since the relative veloc- 
ity of the bodies is increased in binary stars, planetary growth 
will be significantly hi n dered by the pre s ence o f a companion, 
see e.g. iThebault et ail d2006l) ; iThebaultl d2011l) and references 
therein. Here our results indicate that planetesimal growth is 
less negatively influenced because the disk eccentricity is re- 
duced for more realistic radiative disks. As has been shown, ec- 
centric disks tend to increase the mutual relative velocities of 
embedded objects, in particular of differ ent sizes, because of 
the misaligned periastron s of the particles dKlev & Nelsonl2007t 
iPaardekooper et al.l2008l) . Hence, a radiative disk with a low vis- 
cosity could help to promote planetesimal growth. However, it 
remains to be seen how the inclusion of stellar irradiation (from 
both stars) influences the dynamics. Owing to the additional 
heating of the disk, we expect even more mass loss from the 
system and possibly higher disk eccentricities because the disks 
are more isothermal and will have a larger radius. 

Previous studies have indicated that a in mutually in- 
clined system planetesimal growth may be enhanced because 
planetesimals can be size-sorted in differ ently inclined planes 
ilXie & Zhoul 120091; iMarzari et al.l l2009bl) . However, recently 
iFragner et al.l ( 201 ll) showed through full 3D hydrodynamical 
studies that for inclined binaries the relative velocities of dif- 
ferent sized planetesimals increases through inclinations effects. 
They conclude that for inclined systems planetesimal forma- 
tions can take place only for very distant binary stars with 
«bin ^ 60 AU. However, the simulations considered only isother- 
mal disks, and it remains to be seen how radiative effects influ- 
ence the disk. But full 3D radiative simulations are still beyond 
the present computational possibilities, because thousands of or- 
bital timescales of the disk will have to be calculated. 
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